%% 模拟退火算法解决TSP问题
clear,clc                         
%% 初始化
load('city.mat');% 导入城市数据
C=city;% 城市数据两列分别为各个城市的x，y轴坐标
n=size(city,1);% 得到TSP问题的规模,即城市数目
T=100*n;% 设置初始温度
L=100;% 马可夫链长度
K=0.99;% 衰减参数（退火速率）
l=1;% 初始化统计迭代次数为1
len(l)=distance1(city,n);% 每次迭代后的路线长度初始化（目标函数值）
%% 迭代求解 
figure(1); 
e=1e-2;
len3=distance1(city,n);% 迭代过程最优点，不接受劣解，但不影响迭代
while T>e % 停止迭代温度
    len1=len(l);% 计算旧点值
    for i=1:L % 每一个温度下进行L次方向选取
        % 产生新解       
        p1=floor(1+n*rand());% 随机选择两个不同的城市的坐标
        p2=floor(1+n*rand());
        while p1==p2 % 确保选择的两个城市不一样
            p1=floor(1+n*rand());
            p2=floor(1+n*rand());
        end
        tmp_city=city;% 置换两个不同的城市的坐标
        tmp=tmp_city(p1,:);
        tmp_city(p1,:)=tmp_city(p2,:);
        tmp_city(p2,:)=tmp;
        len2=distance1(tmp_city,n);% 计算新路线的距离 
        delta_e=len2-len1;% 新老距离的差值，相当于能量的差
        % 新路线好于旧路线，用新路线代替旧路线
        if delta_e<0        
            city=tmp_city;
            len1=len2;
            if(len2<len3)
                city1=tmp_city;
                len3=len2;
            end
        else
        % 以概率选择是否接受新解
            if exp(-1/T)>rand()% 该处需结合不同的目标函数进行修改，如：
                % exp(-delta_e/T)>rand()
                city=tmp_city;
                len1=len2;
            end
        end
    end
    l=l+1;% 更新迭代次数
    len(l)=len1;% 计算新路线距离 
    T=T*K;% 更新温度
    % 绘制迭代图
    for i=1:n-1
        plot([city(i,1),city(i+1,1)],[city(i,2),city(i+1,2)],'bo-');
        hold on;
    end
    plot([city(n,1),city(1,1)],[city(n,2),city(1,2)],'ro-');
    title(['最短距离:',num2str(len(l))]);
    hold off;
    pause(0.002);
end
figure(2);
for i=1:n-1
    plot([city(i,1),city(i+1,1)],[city(i,2),city(i+1,2)],'bo-');
    hold on;
end
plot([city(n,1),city(1,1)],[city(n,2),city(1,2)],'ro-');
for i=1:n
    text(C(i,1)+0.5,C(i,2),num2str(i));       %标号  
end
text(city(1,1),city(1,2),'    起点');
text(city(n,1),city(n,2),'    终点');
title('最终路线图');
 
figure(3);
plot(len)
xlabel('迭代次数')
ylabel('目标函数值')
title('适应度进化曲线')

% 计算距离的函数，见文件distance1.m
% 从输入City中从上到下依次距离累加，最后计算最后一个点到起点距离，形成一个圈
% function len=distance1(City,n)
% len=0;
% for i=1:n-1
%     len=len+sqrt(sum((City(i,:)-City(i+1,:)).^2));
% end
% len=len+sqrt(sum((City(n,:)-City(1,:)).^2));
% end
% ————————————————
% 版权声明：本文为CSDN博主「Handsome_Zpp」的原创文章，遵循CC 4.0 BY-SA版权协
% 议，转载请附上原文出处链接及本声明。
% 原文链接：https://blog.csdn.net/qq_44865735/article/details/124100778